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Evolution of star clusters in arbitrary tidal fields 
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ABSTRACT 

We present a novel and flexible tensor approach to computing the effect of a time- 
dependent tidal field acting on a stellar system. The tidal forces are recovered from the 
tensor by polynomial interpolation in time. The method has been implemented in a 
direct-summation stellar dynamics integrator (NB0DY6) and test-proved through a set 
of reference calculations: heating, dissolution time and structural evolution of model 
star clusters are all recovered accurately. The tensor method is applicable to arbitrary 
configurations, including the important situation where the background potential is a 
strong function of time. This opens up new perspectives in stellar population studies 
reaching to the formation epoch of the host galaxy or galaxy cluster, as well as for 
star-burst events taking place during the merger of large galaxies. A pilot application 
to a star cluster in the merging galaxies NGC 4038/39 (the Antennae) is presented. 

Key words: globular clusters: general - open clusters and associations: general - 
galaxies: star clusters - methods: analytical - methods: numerical 



1 INTRODUCTION 

The halos of nearly all galaxies are populated by old glob- 
ular clusters that presu mably formed in gaseous discs at 
high redshift (^ ~ 3 - 5, iKravtsov fc GnedinI [20051 '). Young 
'populous clusters', or 'su per-star clust ers' are found in the 
Larg e Magellanic Cloud l|Hode:el Il961h . starburst galaxies 
(e.g. Ivan den Berghl 197ll). interacting galaxie s and merger 



remn ants (e.g. IWhitmore fc Schweizei j 1 19951: 
1997 ) and also in quiescent spirals IjLarsenl 



Miller et al 
2OO2I : iFigei 



20081 ) ■ This suggests that globular cluster formation is 



not unique to the early Universe and that the forma- 
tion of these dense stellar systems is a common phe- 
nomenon in star formation |Elmegreen fc Efremovl 119971 : 



nomenon m star tormation l|±Llmegreen < 
iPortegies Zwart. McMiUan fc Gielesil2010l ) 



To approach this problem numerically is also challeng- 
ing because of the large range of evolutionary time scales 
involved, ranging from several days for tight binaries in the 
cores of clusters to a Hubble time for the h ost galaxy (see 
the recent review bv lDehnen fc Readll201ll ). To be able to 
self-consistently model the evolution of star clusters in 'live' 
gala xies one needs to rely on a direct-tree hybrid approach 
(e.g. iFuiii et al.|[2007l ). 

This is why most studies of the (dynamical) evo- 
lution of star clusters simplify the effect of the exter- 
nal tidal field by a ssuming a static background poten- 
Wein b erg 1990; Vcspcri ni fc Heggiel 
2 00II: iBaumgar dt & M akinol [2003 



tial (e .g. Chernoff 



19971: [Fall fc Z hang 



Dehnen et al|[2bo4,:iGieIes fcBaumgardt„200S: .Hurley et all 



There is increasing evidence from studies of the Milky 
Way and the Andromeda galaxy that some globular clus- 
ters have only recently (past ~ Gyr ) been brought in 



through satellite a ccretion (e.g. JMarm-Franch et al.l l2009l : 
JMackev et al.ll2010r ) . In order to understand the relation be- 
tween the young massive clusters and the old globular clus- 
ters, i.e. their life-cycle, we need to place their evolution in a 
cosmological context. During the formation process of galax- 
ies through repeated accretion phases, substructures such as 
star clusters or dwarf satellites evolve along complex orbits 
in a non-static external potential. This makes their evolution 
difficult (if not impossible) to describe analytically. 

* florent.renaud@cea.fr 



20081: JPeiiarrubia et al.ll2009l : IZonoozi et al.ll201ll ). Ahhough 
this is probably an adequate approximation for many pur- 
poses, it does not suffice for more complicated orbits such as 
those of clusters in mergers of massive galaxies, or satellite 
galajcies that are accreted. 

In light of the last point, several recent studies have 
adopted a (semi-)analytical approach to star cluster evo- 
lution in more realistic extern al tides, such as the h ierar- 
chical build-up of ga laxies (e.g. Prieto fc Gnedinll2008l ) and 
galaxy mergers (e.g. JKruiissen et al.ll2011^ ■ In here the ef- 
fect of mass loss because of stellar evolution, evaporation of 
stars over the tidal radius and the shock-enhanced escape of 
stars because of rapidly varying tidal fields (i.e. disc cross- 
ings and bulge shocks) are applied analytically. In almost all 
cases these processes are assumed to be independent of each 
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other such that the individual resulting mass-loss rates are 
simply added. 

An important, and often dominant, mass-loss mecha- 
nism is the relaxation driven escape of stars, so-called evap- 
oration. Some models assume that a constant fraction of the 
stars evaporate per half-mass relaxation time, independent 
of th e orbit (e.g. iGnedin fc OstrikeJlToOTl : IPrieto fc GnedinI 
120081 ). Others consider that the escape fraction depends on 
the galactoce ntric radiu s , often assumed to be the pericen- 
tre d istance (|King||l962l : llnnanen et"al]|l983l : iFall fc Zhan3 
[2OOII). The final lifetime of clusters can be quite different, 
depending on the details of the assumptions that are made. 

Another critical disruptive agent is mass-loss due to ex- 
ternal tidal perturbations, or 'shocks'. The re lated lifetime 
scales linearly with the density of the cluster (|Spitzeill 19581 : 
lOstriker et al.lll972l ) and it is, therefore, that the results of 
semi-analytical models rely critically on what is assumed for 
the relation between the cluster mass and the half-mass ra- 
dius; a relation which is not only affected by evaporation (i.e. 
because of the reduction of the mas s in time), but also be- 
cause of r elaxation driven expansion (Henon 19651: iGoodmanI 
1 1984 iBaumgardt et al. 2002 ; Gieles ct al. 2010| ), which as 
far as we are aware is neglected in all semi-analytic ap- 
proaches. For clusters on circular orbits in static potentials 
that are well within the tidal limit initially (the half-mass 
radius being less than a few percent of the tidal radius), it 
was found that the expansion phase dominates the evolution 
in the first half of the cluste r's lifetime, while evaporation 
dominates in the second half IjGieles et al.ll201ll ). In the lat- 
ter stage the cluster den sity adjusts to the m ean (galactic) 
density along the orbit (|K{ipper et al.ll2010l ). It is, there- 
fore, necessary to consider both the internal evolution (re- 
laxation) and the external effects (tides) simultaneously if 
one considers the entire life-cycle of clusters. 

Current direct A^-body codes are capable of solving the 
A-body problem numerically for A ~ 10^, together with 
the effects of mass-loss of the individual stars, binary inter- 
actio ns and tidal fields (jPortegies Zwart et al.ll200ll : lAarsethl 
I2OO3I ). Thanks to increasing computational power, it is now 
possible to combine the relaxation driven evolution of star 
clusters in cosmologically motivated external conditions. In 
this paper we present a new method that integrates the ef- 
fect of any tidal field to the evolution of galaxy substruc- 
tures, that we specifically apply to star clusters. We do this 
by extracting the tidal tensor, that contains all information 
about the tidal field at the location of the substructure, from 
galaxy simulations and subsequently 'feed' this to a stellar 
dynamics code. However, we have not yet implemented the 
other half of the scale coupling, as the feedback from the 
small scales (e.g. metal enrichment, stellar winds, super- 
nova explosions) is not retroceded to the ambient galactic 
medium. 

The paper is organised as follows: first, we setup the 
framework for the computation of the tidal acceleration by 
means of tensors (Section [2}. The expressions found are then 
applied to the special cases of circular orbits: well-known 
expressions are retrieved using the new formalism in Sec- 
tion [3l The role of the galactic profile on the evolution of 
clusters through the escape of stars is specially explored 
in Section 13.41 The limitations of the analytical approach 
are explained in Section [l] while Section [5] presents the nu- 
merical implementation of the method to compute the tidal 



forces in a stellar dynamics code. A comparison with pre- 
vious results obtained for idealized configurations is carried 
out. Applications to innovative cases are presented as the 
first practical illustrations of the method. Finally, the lim- 
itations and some possible future developments of our ap- 
proach are discussed. 



2 ANALYTICAL DESCRIPTION OF 
ARBITRARY TIDAL FIELDS 

The main goal of the paper is to provide a general framework 
within which to follow the evolution of self-gravitating stel- 
lar associations in arbitrary and time-dependent tidal fields. 
For concreteness in the remainder of the paper we focus on 
star clusters orbiting within a gala:xy in equilibrium, but the 
formalism can be exported to many other situations (such 
as dwarfs galaxies, galaxy mergers, galaxy clusters, etc). 

2.1 Tidal and effective tensors 

It is convenient to work in coordinates centred on the star 
cluster being embedded in the background gravitational po- 
tential, as opposed to the global system's barycentre. The 
tides derive from gradients in the external gravitational ac- 
celeration across the cluster. Subtracting the acceleration of 
the cluster's centre of mass by the host galaxy, the relative 
acceleration of a member star at the position r' in this frame 
reads 

dV 
df2 



-V<^c(r') - V</.G(r') + V<^g(0), 



(1) 



where (j)^ and (J)q are the gravitational potentials of the 
star cluster and the host galaxy, respectively. The standard 
treatment of the backgro und gravity in the tidal limit (e.g. 
iBinnev fc Tremainell2008h consists in regrouping the last two 
terms in equation dTJ and performing a linear expansion by 
considering that r' <i; 7?g, with _Rg the distance between 
the cluster and the galaxy's barycentre. It is important to 
recall that the linear expansion will hold when and if the 
Laplacian of the galaxy's potential is small at the location 
of the cluster, regardless of the ratio r' /Rq. Bearing this in 
mind, the linearised equations of motion are expressed in 
general form through the tidal tensor Tt of components 



n^ir') 



a'0G 



(2) 



(3) 



(4) 



To first order in r' we have 

V0G(r-') = V<^g(0) - Tt(r') • r' + Oir" 

Substituting in equation ([T]): 

^ = -V<^e(r')+Tt(r')-r'. 

The symmetry T^'' = T^^ allows us to express Tt in diago- 
nal form in the base of its eigenvectors i/^, (i = 1 to 3): the 
amplitude of the eigenvalues A,; is a measure of the strength 
of the tidal field along the corresponding eigenvector. When 
the proper base of the tensor is used to express the accel- 
erations, the reference frame becomes non-inertial. Even so, 
only a rotational component at angular frequency f2 ap- 
pears because the translational component is absorbed in 
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equation ([T|). The net acceleration now includes non-inertial 
terms from fictitious forces: 

gravitational 



dV 

dt2 



internal tidal 

-V(^c(r) +Tt(r)-r 

dn 

Euler 



ft X {ft X r) 



X r — 2n X 



dr 



centrifugal 



_d_f. 
Coriolis 



(5) 



fictitious 

where r is the position vector in the non-inertial frame. The 
centrifugal acceleration can be derived from the gradient of 
a scalar potential 



Mr) = -{r-ilf 



-flV, 
2 



(6) 



defined up to an arbitrary additive constant. This in turn 
leads to an effective tidal potential (l>c and the associated 
effective tidal tensor Te of components 



T:^(r)=r;^(r)+ - 



a'0f 



dx' dx^ 
The total acceleration becomes 



dx^ dx^ 



d^r „ , , ^ ^ , X dn „^ dr 

— = -VUr)+T.{r).r--xr-2nx-. 

In diagonal form, we write the effective tensor Tc as 
Te(r) = 



Ae.l 











Ae,2 











Ac. 3 



(7) 



(8) 



(9) 



with the convention Ao,i ^ Ae,2 ^ Ae,3. In the rest of the 
paper we will refer to the three Ac's as the effective eigen- 
values. 

Equation ^ cannot be simplified for a non-zero Euler 
acceleration. Hence, the following analytical Sections l2.2ll2.3l 
and[3]focus on cases where il is constant in time. More gen- 
eral configurations will be explored numerically in Section [S] 



rt 



/GMc 



1/3 



(10) 



Note that this definition applies to all galactic potentials. 

The sphere of radius rt can be seen as an approximation 
of the physical boundary of the cluster. A more precise three- 
dimensional definition, called the Jacobi surface, can also be 
used. 



2.3 Jacobi surface 

The effective tidal potential derives from the linearization of 
equation (O; 



0c(r) 



r^ ■ Tc(r) ■ r 



(11) 



where r^ is the transpose vector of r. Therefore, with the 
point mass approximation for the cluster potential, the total 
potential is 



{x,y,z) 



GM, Ac,i 

^yx^ + y^ + Z^ ^ 



Ac, 2 2 , Ac, 3 2 

-—y +T—Z 

Ac,l Ac,l 



The three-dimensional surface of equipotential passing 
in Z/i is called the Jacobi surface. From equation (|12|) . we 
find that the corresponding potential energy (also called crit- 
ical energy) is 

E.^^l^. (13) 

2 rt 

The equality of equations p2[) and (|13[) defines the 
equation of the Jacobi surface: 



= 2rt+^/x2 + y^ + z^ ( 



2 , Ac, 2 2 , Ac, 3 2 o 2 1 /, .^ 

X +T-^t/ +-r^2 -3rt).(14) 
Ac,i Ac,i 



.(12) 



A star whose energy is exactly Ej cannot pass through 
this surface, and thus can only escape through the points 
Li or L2, where the surface is 'opened' (see an example in 
Fig. [T]). At energies higher than Ej, the apertures of the 
surface are larger. This plays a non-trivial role in the escape 
rate, as discussed in Section [3.41 



2.2 Tidal radius 

The positions where the internal gravitational acceleration 
of the cluster is exactly balanced by all the other accelera- 
tions are called the Lagrange points Li (with i — 1 to 5) . By 
convention, Li and L2 fall down the galaxy-cluster axis (Li 
being between the two objects). The distance between the 
centre of the cluster and Li is referred to as the tidal radius 
n- 

At Li, it is reasonable to approximate the potential of 
the cluster with that of a point of mass A/c. Furthermore, 
the effective tidal acceleratioiij there is Ae.irt. Finally, the 
Lagrange points are static in this reference frame, meaning 
that the Coriolis acceleration of ii is zero. With these con- 
siderations, equation ^ gives the expression of the tidal 
radius: 



^ The eigenvector related to the largest eigenvalue A(.,i points 
toward the galaxy. 



3 APPLICATION TO CIRCULAR ORBITS 

In this Section, we apply the formulae obtained in Section [2] 
to the special case of circular orbits around various galaxies 
sitting in (— _Rg,0,0). The centrifugal term is derived from 
the rotation speed of the co-rotating reference frame which 
is, by definition, the orbital angular velocity ft — (0, 0, fl): 



fl 



GMgJRg) 
Rl 



(15) 



where Mg[Rg) is the mass of the galaxy enclosed within 
the orbital radius Rg, and Mg (with no argument) is the 
total mass of the galaxy, and M^ <C Mg(-Rg)- 

The centrifugal acceleration is strictly opposed to the 
tidal contribution along the y- and z-axes: 



a'0G 

92/2 



dz2 



(16) 



Rg ^ ^ Rg 

so that the effective tidal eigenvalues are 
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Figure 1. Jacobi surface computed from equation I I14I1 for a clus- 
ter in circular orbit around a point-mass galaxy (equation [Ts)! . 



Ac,i = I - 

Ac,2 = 






Ac. 3 — 



a'0G 



-Rg 



Rg 






Rg 



(17) 



for all circular orbits. 



3.1 Point-mass galaxy 

We first focus on the academic case of a cluster in circular 
orbit around a point-mass galaxy. Using equation HTJI, we 
find that the triplet of effective eigenvalues reads 

GMg 



{Ae,l, Ac, 2, Ac, 3} 



Rh 



■{3,0,-1}. 



(18) 



Note that when replacing this value of Ac.i in equation pop . 
we recover the well-known expression of the tidal radius: 






(19) 



(see iKind Il962l . iFukushige fc Heggid l2000l . or 
iBinnev fc Tremaind I2OO8I ). The corresponding Jacobi 
surface is displayed in Fig. [T] and the one-dimensional 
projections of the total potential (solid black) are compared 
to those of a cluster in isolation (dashed green) in Fig. (2] 
The potential yields a saddle shape in Li (and in L2, by 
symmetry) . 

3.2 Power-law galaxy 

The formalism of the tidal tensor allows us to evaluate the 
tidal acceleration from any galactic potential. As an illus- 
tration, we now focus on power-law galactic profiles of index 
a < 3 whose density is 



-q/2 



(20) 



pG{x,y,z) = po [{x + Rg^ + y^ + z^] 

po being a constant. At the position of the cluster (0,0,0), 
the effective eigenvalues are 



{Ao,i, Ac, 2, Ac, 3} 



47vGpo 



{a, 0,-1}, 



(21) 



(3 — u^-ii-G 

which is comparable to the point-mass case discussed above. 
This sets the value of the tidal radius at 

1/3 



rt = R: 



«/3 



/Mc(3 



/GMc 
Var22 



1/3 



(22) 



The projections of the potential are plotted in Fig. [2] 
for a — 2.5 and 2.0. When normalized to the tidal radius, 
only the z-component differs from the point-mass case. The 
impact of this difference is further explored in Section 13.41 

3.3 Plummer galsixy 

Consider now a IPlummeiJ (|l911l ) potential of characteristic 
radius vq, once again centered on (— _Rg,0,0), 

GMg 



[rf, + (x + Eg)'' + y^ + z^] 



211/2 ■ 



(23) 



We introduce ^ — Rg/to and evaluate the effective eigen- 
values at the position of the cluster: 



{Ao,i, Ac, 2, Ac, 3} = 



GMg 



rui+ef 1(1 +e) 



ie 



,0,-U. (24) 



First, we consider ^ = 1 so that the cluster lies in the 
tidally extensive regime of the Plummer sphere, i.e. where 
the t idal contribution to the first effective eigenvalue is posi- 
tive (jRenaud et al.ll2008l V The triplet of effective eigenvalues 



{Ae,l, Ac, 2, Ac, 3} 



rg \ 8 '^' 4 



which gives a tidal radiu^ of 



Afc 



1/3 



n = ro 



3V2 



J 



Mg 



GM, 



\l^' 



1/3 



(25) 



(26) 



Second, by decreasing ^, we shift the cluster toward the 
centre of the Plummer potential. The tidal contribution to 
Ac,i first tends toward zero; a value reached for ^ = 2~^'2^ 
For smaller values of ^, the cluster lies in the cored region 
of the gal actic potential , and thus is in compressive tidal 
mode (see iRenaud et al.ll2009l ): the tidal acceleration acts 
in the same direction as the intern al gravitational acceler - 
ation of the star cluster. Following IChandrasekharl ([1943), 
Appendix [^demonstrates that the centrifugal contribution 
always compensates the compressive tidal acceleration on 
circular orbits, so that the Lagrange points still exist and 
thus, the tidal radius can be defined, even in compressive 
tidal mode. 

The expression of Ae,i in equation (|24|l reveals that a 
given tidal radius can be obtained at two different galactic 
radii: one in tidally extensive mode, and one in compressive 
mode. The compressive counterpart of our extensive exam- 
ple (equation I26p is obtained for f « 0.66. Both cases are 
plotted in Fig. (2] Interestingly, the potential in the extensive 

2 Substituting a in equation II22I I with the local slope of the den- 
sity profile of the Plummer model (here 5/2 for ^ = 1) does not 
lead to the correct tidal radius because the angular frequency H 
differs between the Plummer and the power-law galaxies. 
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Figure 2. Projections of the total potential (from eQuation ll2l with G = Mc = rt = 1), along the x-, y- and z-axes, for the cluster in 
circular orbit in the x^— plane around a point-mass galaxy, galaxies with power-law density profile of slope 2.5 and 2.0, and a Plummer 
galaxy with the cluster in the tidally extensive (^ = 1) and compressive regimes (^ = 0.66, see text). Along the x- and y-axes, all these 
projections are identical. A circle marks the position of the Lagrange point Li. The potential of the cluster in isolation (r^^, green) is 
also plotted, for comparison. 



case (^ — 1) is strictly identical to that found with a power- 
law galaxy of index a = 1.5: the evolution of identical star 
clusters set in these two configurations would be impossible 
to distinguish. 

The comparison between these five configurations em- 
phasizes that the potential well of a star cluster not only de- 
pends on the tidal radius, but also on the three-dimensional 
shape of the effective tensor, which varies from galaxy to 
galaxy. In the next Subsection, we examine the implication 
of this on the escape rate of stars from the cluster. 



3.4 Escape rate 

For circular or bits, the escape rate varies, to first order, as 



Mc PC n (e.g. iLee fc Ostrikeilll987l : iBaumgardt fc Makinol 



120031 : iGieles fc Baumgardtl 12008). However, the constant of 
proportionality cont ains details of the shape of th e density 
profile of the galaxy. iTanikawa fc Fukushigd (|2010l ) recently 
demonstrated the importance of this secondary effect. They 
did this by considering clusters in galaxies with different 
power-law density profiles (as in our Section [3.2^ . The con- 
stant po was varied such that the tidal radius was the same 
in all cases. They found, somewhat counter-intuitively, that 
in this set-up the clusters with the lowest orbital angular 
velocity Jl had the highest escape rate. In the following, we 
confirm and generalize their result, using the tensor formal- 
ism. 

To escape from the cluster, a star needs (1) to be able 
to fly outside of the Jacobi surface and (2) to exit in a way 
not to fall back in. The first condition implies that the total 
energy E = v'^ /2 + cj> oi the star must exceed the Jacobi en- 
ergy: E > Ej. The second condition tells us that a candidate 
escaper which fulfills the first condition can still be trapped 
in the potential well of the cluster for many crossing-times 
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Figure 3. Projection of the surface of equipotential in the xz- 
plane, for an energy E = 0.95-Ej and for the five galaxies consid- 
ered in Fig. |2] The green circle represent the tidal radius. 



iJFukushige fc Heggid[200ol : lBaumgardtll200ll ): the potential 
barrier keeps increasing with the distance in all directions 
except along the galaxy-cluster axis linking Li and L2 (see 
Fig- El- Therefore, the actual escapers are stars (1) with ex- 
ceeding energy and (2) flying through the apertures in the 
corresponding equipotential surface around Li and L2. The 
size of these apertures depends on the excess of energy and 
on the shape of the tidal field (see Fig. |3]). 
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To evaluate the escape time, we seek the flux of stars 
through the aperture s. We compute it by r epeat ing the cal- 
culation applied by iFukushige fc Heggid (|2000|) to point- 
mass galaxies, but here for any galactic proflle by means of 
the effective eigenvalues. Details of the derivation are given 
in Appendix IbI We find that the timescale for escape for a 
star with an excess of energy E — Ej is 



iesc(-B) 



TvVS 



C GMc 



{E-EjY 



Ae,l 



(27) 



where C is a dimensionless constant which depends on 
the intrinsic properties of the cluster and which we com- 
pute by numerica l integration (C ~ 0.4, see Appendix [Bj. 
iBaumgardd (|200in wrote the (time dependent) dissolution 
timescale tdiss of the entire cluster as 



tdis. oc4/*4/,''(£; = 2B,,), 



(28) 



where i^h is the half-mass relaxation time. This relation 
holds for homologous clusters for which the half-mass radius 
scales linearly with the tidal radius. We can then write the 
dependence of the dissolution timescale for circular orbits 
on the galactic parameters as: 

1/8 



idisi 



3/2 



-^c,3 
Ao,l 



(29) 



We observe a first order dependence of the dissolution 
timescale on the tidal radius, but also a second order effect 
due to the shape of the galactic pote ntial. This relation con- 
firms and extends the conclusion of iTanikawa fc Fukushigd 
lJ2010h to any gala^cy: for a given tidal radius, a highly neg- 
ative ratio of the effective eigenvalues corresponds to a slow 
dissolution, as illustrated in Fig. |4l 

The proportionality factor in equation (|29|) depends on 
the properties of the cluster (i^h, C, Mc), which all evolve 
with time. Thus, tdiss is an instantaneous estimate of the dis- 
solution timescale, and should not be mistaken with the ac- 
tual life-time of the cluster. Moreover, the evolution of these 
properties depends on the tidal field and is very involved 
(if not impossible) to estimate analytically. In Section 15.31 
numerical experiments show indeed that the actual life-time 
can be very different from the analytical value of idisa given 
by equation (|29p . 



4 NON-CIRCULAR ORBITS 

As stated earlier, writing the tidal tensor in its proper base 
makes the reference frame non-inertial: the fictitious accel- 
erations must be added. In the cases of linear or circular 
orbits, the Euler force is null and the formalism of the ef- 
fective tensor allows for a simple mathematical description 
of the net acceleration (see Section [3]). However, for time 
dependent rotations, the Euler effect must be added: the ef- 
fective tensor does not suffice to describe the full impact of 
the galaxy, and the formalism looses its advantage of ana- 
lytical simplicity. 

However, nothing forbids to write the tidal tensor in 
the inertial frame, where it is non-diagonal, and to compute 
the tidal acceleration numerically. The major advantage of 
such a method is that the centrifugal, Euler and Coriolis 
accelerations, which are difficult to evaluate along complex 
orbits, are not required anymore. That is, the tidal tensor 
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Figure 4. Dissolution timescale (tdissi from equation I29|l and 
surface of aperture in tlie Jacobi surface around Li and L2 {A, 
computed from equation I14I I as functions of the ratio of the ef- 
fective eigenvalues, for a given tidal radius and a given cluster. 
Both quantities are normalized to the values they take in the case 
of a point-mass galaxy. The tidally compressive regime is sepa- 
rated from the extensive one by the vertical Une. Symbols mark 
the cases illustrated in Fig. [2] and |3] (from left to right: Plummer 
£, = 0.66, Plummer $ = 1, power-law a = 2.0, power-law a = 2.5, 
point- mass). 



computed in the inertial reference frame fully represents the 
galactic acceleration on star clusters and allows for a numer- 
ical treatment in any galaxy and along any orbit. The next 
Section explains how this method can be implemented in an 
Ai'-body code. 



5 iV-BODY SIMULATIONS: NB0DY6tt 

When expressed in the inertial reference frame, the tidal 
tensor contains all the information on the effect of the galaxy 
on its cluster. Therefore, the equations of motion of all the 
stars of the cluster can be solved numerically, for any tidal 
field. In this Section, we briefly present one implementation 
and a suite of tests, before applying the method to innovative 
cases. 



5.1 Retrieving the external force 

The simulations of the star clusters are done with the stellar 
dynamics code NB0DY6 and its version for Graphics Process- 
ing Unit (GPU, Aarseth 201C|j). Several cases of galactic po- 
tentials already exist in NB0DY6 and have been widely used 
in previous studies. For testing purposes, we use these fea- 
tures and refer to them as built-in methods. On top of these 
pre-existing tools, we have modified NBDDY6 to include the 
tidal forces by means of the tidal tensor: this new version of 
the code is called NBDDY6tt. 

Before the simulation, the tidal tensor is computed in 
the inertial reference frame (equation [2]): its nine compo- 
nents are sampled along the orbit of the cluster within the 



|http : //www . ast . cam . ac ■ uk/~sverre/web/pages/nbody ■ htm | 
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galaxy, either analy tically, or by means of independent galac- 
tic simulations (see lRenaud et al.ll2009l : lRenaudll201Cl ). The 
sampling frequency is chosen to ensure that the high fre- 
quency features of the tidal acceleration, both in term of 
intensity and orientation, are recovered. A table of sampled 
tensors is then passed to NB0DY6tt. During the simulation of 
the cluster, the components of the tensor are quadratically 
interpolated whenever the gravitational force on a particle 
needs to be updated. For a star of mass m at the position 
{a;'*} with respect to the centre of the cluster, the tidal force, 
whose i-th component reads 



K 



mY^Tl^x'^ 



(30) 



is added to t he gravitation al force due to the A'' — 1 other 
particles (see I Aarsethll2003l , for details on the solving of the 
A-body problem). 

As soon as physically-time-dependent processes (like 
stellar evolution) are not involved, the entire study remains 
scale-free. The units adopted below are the cluster's A-body 
units (G — Mc = — 4_Ec = 1, with Ec being the total energy 
of the cluster; see lHeggie fc MathieulligSQ l. One of the pos- 
sible scalings to physical units is proposed in Table [T] 

In the inertial reference frame, it is generally not possi- 
ble to evaluate the effective eigenvalues and thus, the tidal 
radius and the energy of the Jacobi surface. Therefore, to 
define the cluster membership, we have adopted an empiri- 
cal criterion; a star is considered as a cluster member when 
the sum of its kinetic energy exactly balances the potential 
energy due to the A — 1 other s tars, A being determin ed it- 
eratively until convergence f see IPenarrubia et aLlbOOGf ). For 
circular orbits, we have measured that this count only differs 
from the number of stars within the tidal radius by a few 
percent. 
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Figure 6. Number of stars in tlie cluster along the orbits A, B 
and C around a point-mass galaxy, normalized to its initial value. 
The black curves represent the solution using the built-in method 
of NBDDY6, while the coloured ones are associated with the use of 
the tidal tensor. The symbols mark the passage of the cluster at 
the positions marked in the top-right panel of Fig. \5\ 



of the two approaches is very good for all the orbits. In 
particular along the elliptical orbit C, the expansion of the 
outermost layers of the cluster and the increased mass-loss 
near the pericentre passages is well-reproduced by the new 
method, both in term of time (epoch, delay and duration) 
and amplitude. These tests demonstrate that the interpola- 
tion scheme used to evaluate the tidal tensor at any time 
and the computation of the force done by NBDDY6tt allow 
us to retrieve the results obtained with well-tested methods, 
at a high level of accuracy. 



5.2 Tests 

To ensure the validity of the method, we ran a series of 
tests to compare the new implementation with the built- 
in methodcl of NB0DY6. For both codes, the centre of the 
cluster is fixed and the tidal field mimics the orbit of the 
galaxy around it. The orbits, shown in the top-right panel 
of Fig. O are circular of radius Ro — 1000 (A) , circular of 
radius Rg = 3000 (B), and elliptical with an eccentricity 0.5 
and a pericentre of 1000 (C), which places the apocentre at 
a distance of 3000 so that the tidal field takes intermediate 
values, between the s trong field (A) a nd the weaker one (B). 
The cluster is a IPlummerl l|l91H ) sphere made of A = 
8000 equal-mass particles. To focus on the effect of the tides, 
we have switched off stellar evolution. The parameters of the 
run are listed in Table [l] The evolution of some Lagrange 
radii of the cluster are plotted in Fig. [5] before core-collapse, 
the relative differences (1 — ''NB0DY6/^NB0DY6tt) remains be- 
low 5%, in all cases. The differences increases after core col- 
lapse, but not systematically, probably because of the forma- 
tion of binaries which is sensible to numerical and A-body 
noises. As a complement. Fig. |6] plots the evolution of the 
number A of stars in the cluster, normalized to its initial 
value. In both measurements ('"Lagj. and A), the agreement 



Option #14 = 3 in NBDDY6. 



5.3 Other galactic profiles 

The method having been successfully tested, we now ex- 
plore the evolution of clusters, still on circular orbits, but 
in the galactic potentials presented in Section [3] and Fig. (2] 
The parameters Mq and po have been chosen so that the 
tidal radius is the same as that of the orbit A. The ten- 
sors are computed analytically along the orbit and passed 
to NB0DY6tt. 

The evolution of A is presented in Fig. [T] The life-times 
of the clusters are ordered as predicted in Section 13.41 We 
notice however that, for example, the life-time of the cluster 
in the Plummer galaxy at ^ = 0.66 is « 2.2 time longer than 
that around the point-mass galaxy, while equation (|29|) pre- 
dicts a value of ~ 1.57. The discrepancy is due to a different 
evolution of the internal properties of the cluster (trh, C, 
Mc), which are enclosed in the constant of proportionality 
of equation (|29p . Furthermore, because of the different shape 
of the potential (i.e. the second order effect in equation I29|l . 
the escape of stars does not occur at the same rate in all 
galaxies, which modifies Ale and the tidal radius itself (i.e. 
the first order effect) differently from galaxy to galaxy. As 
a consequence, the effect due to the mild variation of the 
three-dimensional shape of the Jacobi surface gets strongly 
amplified with time. 

Simulations at higher resolution (A = 16 x 10'' and 32 x 
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Table 1 . One possible scaling of the simulations in physical units 



Quantity 


A''-body units 




Physical units 


(at initial time) 










Cluster scale 






Mass (Mc) 


1 




8 X 10^ M0 


Mass of a particle (m) 


1/8000 




IMq 


Virial radius (r^) 


1 




Ipc 


Characteristic radius (ro,c) 


37r/16 




0.59 pc 


Half-mass radius (rj,) 


ro,c/V22/3 - 1 




0.77 pc 


Crossing time (tcr) 


2V2 




0.47 Myr 




Galactic scale 






Mass (Mq) 


1.25 X 10^ 




lO^O Mq 


Y'^ effective eigenvalue* (Ao,i) 


[3.75,0.13] X 10-3 


[135.C 


,5.0] X 10-3 Myr-2 


Orbital radius* {Rq) 


[1, 3] X 10^ 




[1,3] kpc 


Orbital pcriodt {tarh) 


[178.3,927.7,504.8] 


[29 


6,154.0,83.8] Myr 



* for the orbits A and B, respectively. 
t for the orbits A, B and C, respectively. 
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Figure 7. Evolution of the normalized number of stars in the 
cluster on circular orbit in the galactic potentials presented in 
Section [3] compared to the case of a cluster in isolation (dashed 
green line). 



10 particles, but always the same Mc) have also been done, 
with these tidal fields. Not surprisingly, NB0DY6tt shows that 
the actual life-time of the clusters increases as ~ A'^'^' * (see 
equation 1281 with the usual scaling of trh with A''), so that 
the discrepancy between the analytical values of tdiss and the 
life-times measured numerically for our five galactic cases is 
preserved at higher resolution. This demonstrates that our 
conclusions are not affected by the low- A*' statistics of our 
experiments. 

It is therefore very involved to derive analytically the 
mass or the life-time of clusters from initial parameters only. 
The use of numerical methods like NB0DY6tt provides a so- 
lution, but at the non-negligible cost of computational time. 



5.4 Fully arbitrary tidal field 

As a last step toward generality, this Subsection presents the 
results obtained for a complex and highly time-dependent 



tidal field. The orbit chosen is extracted from a simulation 
of the Antennae galaxies (NGC 4038/39), a prototypical ma- 
jor merger. The galactic run, the para meters and the orbit 
are described in iRenaud et al.l (120091 . Fig. 8, orbit B, see 
also their Fig. 3): the cluster starts orbiting in the disc of 
NGC 4038 at ~ 6 kpc (on average) from the galactic cen- 
tre; then it is ejected by the first galactic pericentre pas- 
sage into the intergalactic bridges before falling back into 
the central region and remaining there for the rest of the 
merger. The Antennae being a real, observed object, the 
units are now scaled according to the galactic simulation 
(based on the observed spatial extension of the tidal tails 
and the peak radial velocity). The NB0DY6tt run is arbitrar- 
ily started 100 Myr before the first pericentre passage of the 
two galaxies. The orbit of the cluster has also been inte- 
grated within its host galaxy (NGC 4038) in isolation, as a 
reference simulation. In both cases, the cluster is setup iden- 
tically to those of the previous Sections (see Table [TJ . Our 
physical scaling makes it comparable in mass (8000 Mq) 
and density (~ 2000 Mq pc"*^ within the half-mass radius ) 
to Westerlund 1 or NGC 3603 (|Portegies Zwart et al-lboiOl ). 
The maximum eigenvalue Ai of the tidal tensor is plot- 
ted in Fig. |8]a, in the case of the merger (red) and of the 
isolated galaxy (blue). The cluster is in compressive mode 
for Ai < 0. The centrifugal terrr|j f2^ is shown in logarith- 
mic scale in Fig.[8]b. The peaks denote the velocity kicks the 
cluster receives when it is gravitationally slingshot. Fig. [Sjc 
displays the ratio of the density of the cluster and the local 
density of the galaxjij. Finally, the evolution of the number 
of stars in the cluster is displayed in Fig. |Sld, and compared 
to that of the cluster in isolation (green). 

^ To evaluate the centrifugal term O^, we use the eigenvector vi 
associated with Ai. This vector points towards the main source 
of gravitation and thus, the variation of its direction gives the 
instantaneous orbital rotation speed. We obtain it via the dot 
product of two consecutive (normalized) i/i's: 



n(t) : 



acos[i'l(t) • i'i(t — dt)] 
dt ■ 



° The local density of the galaxy is given by the trace of the tidal 
tensor, through Poisson's law: 
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Figure 5. Top-right panel: orbits of the clusters in the galactocentric frame. Other panels: Evolution of the 10%, 50% and 90% Lagrange 
radii of the cluster on the three orbits, computed by the built-in method (black) and by NB0DY6tt (colour). On the top of each panel, 
crosses mark the passages of the cluster at its left-most position in the top-right panel (i.e. at pericentre for orbit C). The core-collapse 
phase is well visible at t ~ 1700 - 2000. 



At i ~ 180 Myr, the cluster is in the extensive regime of 
the tidal bridges, marked by a slow increase of Ai: the mass- 
loss accelerates with respect to that of the cluster in the 
isolated galaxy. During the second galactic pericentre pas- 
sage (i ~ 285 Myr), the tidal field is mostly extensive, which 
once again, enhances the escape of stars. Later, the galaxies 
have merged and the cluster is orbiting the remnant. Its or- 
bital period (~ 140 Myr) is clearly visible in both Ai and Q,^, 
but also in the mass-loss which is accelerated at pericentre, 
in a comparable fashion as the elliptical test-case of Fig. |6l 
The rapid variations of Ai near the pericenter passages exist 
because the cluster flies at high speed in a highly asymmet- 
ric potential, which can be momentarily compressive at the 
location of the cluster. 

Along this particular orbit, the local galactic density is, 
on average, about three orders of magnitude smaller than 
that of the cluster but the ratio of the two is peaked at sev- 
eral epochs. Therefore, the tidal field only affects the clus- 
ter during precise and short periods of time: the cluster is 
dense, robust enough to remain mildly disturbed on the long 



E- 



-\''^(j,G(r) = -4nGpG(r}. 



timescale. In other words, the life-time of this cluster is al- 
most independent of the merger; only its 'life-style' differs 
from the case of the isolated galaxy. However, the analysis of 
this single case is not statistically relevant to conclude that 
the merger has no secular impact on its clusters. The prop- 
erties of the cluster population, in particular the cluster age 
function, rely on many parameters (structural and orbital) 
and we leave their study to a forthcoming paper. 



6 SUMMARY, LIMITATIONS AND 
CONCLUSIONS 

In this contribution, we propose a new formalism to describe 
the tidal field in A'^-body simulations by the means of ten- 
sors. Although the analytical approach rapidly becomes too 
involved for accelerated motions, we have derived the ex- 
pressions of quantities representing the effect of the tides 
on stellar systems for circular orbits, with no restriction on 
the shape of the external potential. The main results of this 
study are: 

• The use of the tidal and effective tensors allows us to 
simplify the representation of the problem, without loss of 
information (Section l2.ip . Useful quantities like the tidal ra- 
dius (equation llOp , the energy of the Jacobi surface (equa- 
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Figure 8. (a): Maximum eigenvalue of the tidal tensor along the 
orbit in the Antennae (see text, red), and in the same progenitor 
galaxy but isolated (blue), (b): Instantaneous centrifugal contri- 
bution to the effective acceleration, (c): Logarithm of the ratio 
of the density of the cluster and the local density of the galaxy. 
The horizontal dashed lines show the ratios averaged over the 
life-time of the cluster, (d): Evolution of the normalized number 
of stars in the clusters in these galaxies, and in an isolated cluster 
(green). Vertical dotted lines mark the galactic pericentre pas- 
sages and the arrow indicates when the cluster is in the bridges 
in the simulation of the Antennae. 



tion ll3|l or the escape timescale (equation I27|) can be easily 
computed for a given cluster. 

• The tidal radius is a first order approximation of the 
effect of a galaxy on a star cluster. The three-dimensional 
shape of the galactic potential has a second order effect (Sec- 
tion [213] and equation I29p . 

• For a given tidal radius, a point-mass is the most ef- 
ficient galactic profile in dissolving a cluster. Shallower po- 
tentials lead to longer dissolution times (Fig. Q. 

• A cluster in compressive tidal mode loses its stars more 
rapidly than in isolation, but significantly more slowly than 
if it was in extensive mode (Section 13. 4p . 



• The knowledge of the evolution of the cluster parame- 
ters (half-mass radius, mass, energy distribution) is key to 
estimate the cluster life-time (Section l3.4l and Appendix IB)) . 

• For non-circular orbits, the analytical approach be- 
comes very involved so that general and simple formulae 
do not exist (Section [Jjl. 

To overcome this issue, we have developed and implemented 
a numerical method called NBDDY6tt which computes the 
evolution of A'^-body models of star clusters in any tidal 
field, by means of the pre-calculation of inertial tidal ten- 
sors. This method has been successfully tested and applied 
to innovative cases. Our main conclusions are: 

• The tidal force felt by the stars of a star cluster can 
be accurately computed by evaluating the tidal tensor at all 
time by means of quadratic interpolation (Section 15. 2p . 

• The dependence of the dissolution time on the three- 
dimensional shape of the galactic potential highlighted in 
Section [3.41 is confirmed by our numerical experiments. The 
numerical method does not suffer from the lack of knowledge 
of the evolution of the cluster properties that one has to face 
in a (semi-) analytical approach (Section 15. 3|l . 

• Our implementation has been applied to the complex 
case of a cluster in the major merger of the Antennae galax- 
ies (Section l5.4|l . Its mass-loss refiects the nature of the time- 
dependent tidal field experienced along the orbit. In particu- 
lar, the alternation of extensive and compressive tidal modes 
strongly affects the instantaneous dissolution rate, over a 
time-scale of several 10^ yr. 

This preliminary study shows the way to a very wide 
range of possible applications. However, one should keep in 
mind that our study is limited in several aspects. On the one 
hand, our cluster simulations are gas-free. Therefore, we can- 
not address the important point of the early life of the star 
cluster, prior to gas expulsion (first ~ 10^ yr of the cluster's 
life), and we limit our study to initially relaxed systems. On 
the other hand, the second half of the galaxy-cluster cou- 
pling is not taken into account by our method: the feedback 
from stellar evolution, the escape of stars in the field and the 
formation of long tidal tails or streams are not implemented. 
Although they have a limited impact on the evolution of the 
cluster itself, it would be important to monitor their effects 
at larger scale (e.g. the metal enrichment of the interstel- 
lar medium). Furthermore, the mass- loss experienced by a 
cluster could affect its orbit within its host galaxy and thus, 
change the tidal field. In our approach, the scale decoupling 
does not allow to account for such effect. 

To conclude, we have shown that a star cluster plunged 
in a time-dependent tidal field leads to a complex evolu- 
tion, out-of-reach of (semi-) analytical approaches. A numer- 
ical method like the one proposed by NB0DY6tt provides a 
framework for future explorations of the role of the tides 
on star clusters. Among others, the use of a stellar mass 
function, primordial binaries and high order stars, stellar 
evolution and stellar mass-loss are as many lines of investi- 
gation to be pursued on top of the evolution in a background, 
time-dependent, galactic potential. 

In a forthcoming paper, we will use more detailed galac- 
tic simulations including a prescription on the star forma- 
tion, to compute the tidal tensors of an entire cluster popu- 
lation. NB0DY6tt will help us to derive the cluster- mass and 



Evolution of star dusters in arbitrary tidal fields 11 



age functions and their evolution, in various type of galaxies, 
in particular in mergers. 
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APPENDIX A: EXISTENCE OF THE 
LAGRANGE POINTS 

The Lagrange points Li and L2, which define the tidal radius 
rt, exist when the internal gravitational acceleration of the 
cluster is balanced by the effective tidal acceleration: 

GMc 



Ac, in. 



Along a circular orbit of radius Rq , one gets 



GMc 






+ n" 



rt, 



(Al) 



(A2) 



which can be re-written by introducing the epicycle fre- 
quency k: 

GMc 



= [- (k^ - 30^) + n^] 



rt- 



(A3) 



This tells us that the tidal radius can be defined for 
K^/Q,'^ < 4, which is always the case since the maximum 
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value K^ /Q? = 4 is reached for homogeneous mass distribu- 
tions. To conclude, on a circular orbit, the centrifugal ac- 
celeration always compensates the tidal component, even in 
compressive mode, so that the Lagrange points Li and I/2 
always exist. 



APPENDIX B: COMPUTATION OF THE 
ESCAPE TIME 

In this Appendi x, we generalize the e x press ion of the escape 
time derived bv iFukushige fc Heggid (|200(1 ) for a cluster in 
a circular orbit around a point-mass galaxy, to the case of 
any galactic potential. First, the origin of the coordinates is 
shifted to L\ and the total potential is expanded to second 
order, so that 



<}>-Ej 



Ae,l 

~2~ 



-3a:: +y + z 1 



An -1 



(Bl) 



The flux of phase volume across the new x = is expressed 
as 



HE) = / 5 



I -|- — E ) X Ax Aij Az Ay Az , 



(B2) 



where the dot indicates derivation with respect to time and 
5 is the Dirac function. We change variables to ui : i 1— >■ 
(j) -\- v^ /2 ~ E so that Aw = xAx and integrate the Dirac 
function over w to get 



T{E) = / AyAiAyAz, 
J ±>o 

with integration boundaries satisfying 



2{E - E3) - Ac,i 



y^ + z^\l 



Ac 



>0. 



(B3) 
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We change the remaining four variables into the hyper- 
spherical coordinates {i?, S, r, t/j}, i.e. 



, -1/2 n 
y = A^ / _Rcos 



Ae,l(l-^ 



y — R sin 9 sin r cos ■;/' 
z = R sin sin r sin ip 



-1/2 



7? sin 6 cos r 



with ■ 



R>0 

e € [0, 27r] 
r G [0, -k] 
V' G [0, tt] 



so that the condition equation (|B4[I becomes 

2{E - Ej) -R^ >0. (B6) 

The determinant of the Jacobian matrix of the transforma- 
tion gives the hyper- volume element: 

R^ sin^ 9 sin r 



dydidyd^; 
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The flux is finally 

2n' {E-Ejf 
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The total flux is 2J-" because stars can escape through aper- 
tures around two Lagrange points. 

Similarly, the phase-space volume can be written 



V 



<5U + 



E I A^rA^v. 



(B9) 



We first take out the angular part of the velocity and inte- 
grate the Dirac function over v. 



V = An y/2(E~(t>)d^r. 



(BIO) 

(Bll) 

(B12) 
(B13) 



is a dimensionless quantity which describes the intrinsic 
properties of the cluster. Instead of integrating over the en- 
tire solid angle, we note that the fiux V is non-zero only 
at the vicinity of the Lagrange points Li and L2 . In a first 
approximation, we may consider a non-zero flux only at the 
exact position of Li and L2. In that case, we replace the an- 
gular dependence of the previous integral with Dirac func- 
tions so that only two angular directions remain. That is, C 
becomes 



Defining the dimensionless quantities 

** ={E- (j>) rt/iGMc) 
r* — r/rt, 

and substituting in equation (|B10|) yields 

V = 47r^ (GMc)^''^ r"^^^ C, 

where 



C= / %/^dV* 



C = 2 



/** 



■ Ar* 



(B14) 



For King profiles with the usual parameter *I'o/(t^ ranging 
from 3 to 12, we found (by means of numerical integrations) 
that C takes values ranging from 0.38 to 0.39 with a maxi- 
mum reached for ^o/o"^ ~ 8. 

By using the general definitions given in Section [2j we 
find 

V = 47r CV2 {GM,)*/''X;l^\ (B15) 

It follows that the timescale for escape with the energy E is 

V 



t.sc{E)- 

(B5) 
which gives equation (|27 



(B16) 



